---
title: "Figures for Mapping violence perceptions through YouTube comments"
format: html
editor: visual
---

## Figure 1

```{r}
library(tidyverse)
library(slider)
library(sf)
library(cowplot)
library(terra)

#### Data prep ####

mx_municipal_crime_stats_2025 <- 
  read.csv("data/mx-municipal-crime-stats-2025.csv", encoding = "latin1") |>
  dplyr::filter(`Subtipo.de.delito` == 'Homicidio doloso') |>
  tidyr::pivot_longer(
    cols = Enero:Diciembre,  # Select all month columns
    names_to = "mes",
    values_to = "value"
  ) %>%
  dplyr::mutate(
    # Convert Spanish month names to numbers
    month_num = case_when(
      mes == "Enero" ~ "01",
      mes == "Febrero" ~ "02",
      mes == "Marzo" ~ "03",
      mes == "Abril" ~ "04",
      mes == "Mayo" ~ "05",
      mes == "Junio" ~ "06",
      mes == "Julio" ~ "07",
      mes == "Agosto" ~ "08",
      mes == "Septiembre" ~ "09",
      mes == "Octubre" ~ "10",
      mes == "Noviembre" ~ "11",
      mes == "Diciembre" ~ "12"
    ),
    # Create YYYY-MM column
    month = paste(Año, month_num, sep = "-")
  ) %>%
  dplyr::select(-mes, -month_num) |>
  dplyr::group_by(`Cve..Municipio`, month) |>
  dplyr::summarise(value = sum(value))
  
mx_tot <- 
  mx_municipal_crime_stats_2025 |>
  dplyr::group_by(month) |>
  dplyr::summarise(tot = sum(value)) |>
  dplyr::mutate(month = 
                  as.Date(paste(month, "01", sep = "-")),
                tot_ma3 = slide_dbl(tot, mean, .before = 1, .after = 1, .complete = TRUE),
                tot_ma5 = slide_dbl(tot, mean, .before = 2, .after = 2, .complete = TRUE))

mx_2020_24_mean <-
  mx_municipal_crime_stats_2025 |>
  dplyr::mutate(month = as.Date(paste(month, "01", sep = "-"))) |>
  dplyr::filter(month >= as.Date("2020-01-01") & 
                  month < as.Date("2025-01-01")) |>
  mutate(
    # Ensure Cve..Municipio is numeric
    Cve..Municipio = as.numeric(Cve..Municipio),
    # Format to 5-digit code with leading zeros, then add MX prefix
    ADM2_PCODE = paste0("MX", sprintf("%05d", Cve..Municipio))
  ) |>
  dplyr::group_by(ADM2_PCODE) |>
  dplyr::summarise(mean = mean(value))

mx_2024_diff <- 
  mx_municipal_crime_stats_2025 %>%
  filter(month %in% c("2023-05", 
                      "2024-05")) %>%
  mutate(
    # Ensure Cve..Municipio is numeric
    Cve..Municipio = as.numeric(Cve..Municipio),
    # Format to 5-digit code with leading zeros, then add MX prefix
    ADM2_PCODE = paste0("MX", sprintf("%05d", Cve..Municipio))
  ) |>
  select(ADM2_PCODE, month, value) |>
  pivot_wider(names_from = month, 
              values_from = value,
              names_prefix = "murders_") %>%
  mutate(
    murders_diff = `murders_2024-05` / `murders_2023-05`,
    murders_2023 = `murders_2023-05`,
    murders_2024 = `murders_2024-05`
  ) %>%
  select(ADM2_PCODE, murders_2023, murders_2024, murders_diff)
             

mx.sf <- 
  read_sf("data/mx-municipal-shape/mex_admbnda_adm2_govmex_20210618.shp")

mx_pop <- 
  read.csv("data/mex_admpop_adm2_2024.csv", encoding = "latin1")


mx_2020_24_mean$pop <- 
  mx_pop$T_TL[match(mx_2020_24_mean$ADM2_PCODE,
                    mx_pop$ADM2_PCODE)]

mx.sf$mean <- 
  (mx_2020_24_mean$mean /
  (mx_2020_24_mean$pop / 100000))[
    match(mx.sf$ADM2_PCODE,
          mx_2020_24_mean$ADM2_PCODE)
  ]

mx.sf$diff <- 
  mx_2024_diff$murders_diff[
       match(mx.sf$ADM2_PCODE,
             mx_2024_diff$ADM2_PCODE)
     ]


gg_map_2020_24 <- 
  mutate(mx.sf, 
         value_category = cut(mean, 
                              breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.5, 5, Inf),
                              labels = c("0-0.2", "0.2-0.4", "0.4-0.6", "0.6-0.8", 
                                         "0.8-1.0", "1.0-1.5", "1.5-2.5", "2.5-5", "Over 5"),
                              include.lowest = TRUE,
                              right = FALSE)
  ) |>
  ggplot() + 
  geom_sf(aes(fill = value_category), color = NA) +
  scale_fill_brewer(palette = "YlOrRd", direction = 1) +
  labs(
    title = "Monthly Average Homicide Rate by Municipality per 100,000 people",
    subtitle = "2020-2024 (incidents per month)",
    fill = NULL,
  ) +
  theme_void() +
  guides(fill = guide_legend(nrow = 1)) +
  theme(legend.position = "bottom")

gg_timeline <- 
  ggplot(data = mx_tot) +
  geom_line(aes(x = month, y = tot_ma5, color = "5-month moving average", linetype = "5-month moving average")) +
  geom_line(aes(x = month, y = tot,  color = "Monthly average", linetype = "Monthly average"), alpha = .5) +
  geom_point(aes(x = month, y = tot, color = "Monthly average"), alpha = .2) +
  geom_vline(xintercept = as.Date("2024-06-02"), 
             color = "red", 
             linetype = "dashed", 
             size = 1) +
  annotate("text", 
           x = as.Date("2024-06-02"), 
           y = max(mx_tot$tot, na.rm = TRUE) * 1.01,
           label = "General Election\n(2 June 2024)", 
           hjust = -0.1, 
           color = "red",
           size = 3) +
  labs(
    title = "Monthly Nation-wide Timeline with 5-Month Centered Moving Average",
    x = NULL,
    y = "incidents",
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_color_manual(
    name = NULL,
    values = c("5-month moving average" = "black", "Monthly average" = "red")
  ) +
  scale_linetype_manual(
    name = NULL,
    values = c("5-month moving average" = "solid", "Monthly average" = "solid")
  ) +
  theme(
    legend.position = "bottom",
    legend.key = element_rect(fill = "white")
  )

# Figure 1

ggsave(filename = "figures/fig-1.jpg", width = 12, height = 12,
       plot_grid(gg_map_2020_24, 
                 gg_timeline,
                 nrow = 2)
)

ggsave(filename = "figures/fig-1.pdf", width = 12, height = 12,
       plot_grid(gg_map_2020_24, 
                 gg_timeline,
                 nrow = 2)
)

```

## Figure 3

```{r}
library(tidyverse)
library(cowplot)
library(scales)

load("data/fig-3.RData")

cowplot_by_day_ts <- 
  cowplot::plot_grid(
    
    ggplot(videos_by_day, 
           aes(x = date, y = n)) +
      geom_line() +
      scale_x_date(
        date_labels = "%b %Y",           # Format: Jan 2020
        date_breaks = "3 months",        # Break every 3 months
        limits = c(as.Date("2020-01-01"), ("2024-06-01"))) +
      labs(y = "n. videos", x = NULL) +
    
        geom_vline(xintercept = as.Date("2020-03-10"), linetype = "dashed", color = "darkgrey", alpha = 0.7) +
  geom_vline(xintercept = as.Date("2020-06-05"), linetype = "dashed", color = "darkgrey", alpha = 0.9) +
  geom_vline(xintercept = as.Date("2020-07-15"), linetype = "dashed", color = "darkgrey", alpha = 0.9) +
  geom_vline(xintercept = as.Date("2023-01-05"), linetype = "dashed", color = "darkgrey", alpha = 0.9) +
  geom_vline(xintercept = as.Date("2024-06-02"), linetype = "dashed", color = "darkgrey", alpha = 0.9) +
  # Event annotations with segments pointing to dates
  # March 2020 - Women's protests against violence (MOVED ABOVE SPIKE)
  annotate("text", 
           x = as.Date("2020-03-10"), y = 10500,
           label = "Mar 2020\nMassive protests\nagainst violence", 
           size = 2.8, hjust = 0.5) +
  # June 2020 - Guadalajara protests
  annotate("text", 
           x = as.Date("2020-06-05"), y = 8000,
           label = "Jun 2020\nGuadalajara\nprotests", 
           size = 2.8, hjust = 0.5) +
  # July 2020 - Guanajuato crime wave
  annotate("text", 
           x = as.Date("2020-07-15"), y = 4000,
           label = "Jul 2020\nGuanajuato\ncrime wave", 
           size = 2.8, hjust = 0.5) +
  # January 2023 - Sinaloa unrest
  annotate("text", 
           x = as.Date("2023-01-05"), y = 7500,
           label = "Jan 2023\nSinaloa unrest\n(Guzmán arrest)", 
           size = 4, hjust = 0.5) +
  # 2024 Election campaign box
  annotate("text", 
           x = as.Date("2024-03-15"), y = 1000,
           label = "2024 Federal election", 
           size = 4, hjust = 0.5)
      
      ,
    
    ggplot(comments_by_day, 
           aes(x = date, y = n)) +
      geom_line() +
      scale_x_date(
        date_labels = "%b %Y",           # Format: Jan 2020
        date_breaks = "3 months",        # Break every 3 months
        limits = c(as.Date("2020-01-01"), ("2024-06-01"))) + 
      labs(y = "n. comments", x = NULL) +
    
            geom_vline(xintercept = as.Date("2020-03-10"), linetype = "dashed", color = "darkgrey", alpha = 0.7) +
  geom_vline(xintercept = as.Date("2020-06-05"), linetype = "dashed", color = "darkgrey", alpha = 0.9) +
  geom_vline(xintercept = as.Date("2020-07-15"), linetype = "dashed", color = "darkgrey", alpha = 0.9) +
  geom_vline(xintercept = as.Date("2023-01-05"), linetype = "dashed", color = "darkgrey", alpha = 0.9) +
  geom_vline(xintercept = as.Date("2024-06-02"), linetype = "dashed", color = "darkgrey", alpha = 0.9)
    ,
    
    nrow = 2
  )


cowplot_by_24h <-
  cowplot::plot_grid(
    
    ggplot(videos_by_24h, aes(x = time_posix, y = n)) + geom_line() +
      scale_y_log10() +
      scale_x_datetime(
        date_labels = "%H:%M",
        date_breaks = "2 hours"
      ) +
      annotation_logticks() +
      labs(x = NULL, y = "n. videos"),
    
    ggplot(comments_by_24h, aes(x = time_posix, y = n)) + geom_line() +
      scale_y_log10() +
      scale_x_datetime(
        date_labels = "%H:%M",
        date_breaks = "2 hours"
      ) +
      annotation_logticks() +
      labs(x = "Mexico City time", y = "n. comments"),
    
    
    nrow = 2
  )


cowplot_grid_pop <- 
  cowplot::plot_grid(
    
    ggplot(r_df, aes(x = value, 
                     y = n_videos)) +
      geom_point(alpha = .5, size = .8) +
      scale_x_log10(limits = c(100, NA),
                    labels = label_comma()) +
      # scale_y_log10(labels = label_comma()) +
      annotation_logticks() +
      labs(x = "population (log.)", y = "n. videos"),
    
    ggplot(r_df, aes(x = value, 
                     y = n_comments)) +
      geom_point(alpha = .5, size = .8) +
      scale_x_log10(limits = c(100, NA),
                    labels = label_comma()) + 
      scale_y_log10(labels = label_comma()) +
      annotation_logticks() +
      labs(x = "population", y = "n. comments"),
    
    nrow = 2
  )


bottom_row <- 
  cowplot::plot_grid(cowplot_by_24h, 
                     cowplot_grid_pop,
                     ncol = 2)


ggsave(filename = "figures/fig-3.jpg", width = 14, height = 10,
       
       cowplot::plot_grid(cowplot_by_day_ts, 
                          bottom_row,
                          nrow = 2,
                          rel_heights = c(.7, 1))
       
)

ggsave(filename = "figures/fig-3.pdf", width = 14, height = 10,
       
       cowplot::plot_grid(cowplot_by_day_ts, 
                          bottom_row,
                          nrow = 2,
                          rel_heights = c(.7, 1))
)
```

## Figure 4

```{r}
library(tidyverse)
library(cowplot)

load("data/fig-4.RData")

fig_4 <- 
  plot_grid(
         
         ggplot(r_df, aes(x=x, y=y)) + 
           geom_tile(aes(fill = clean_labels)) +
           scale_fill_brewer(palette = "Spectral", direction = 1) +
           theme_void() +
           labs(
             title = "Population 2020",
             fill = NULL,
           ) +
           guides(fill = guide_legend(nrow = 2)) +
           theme(legend.position = "bottom",
                 plot.background = element_rect(fill = "white", color = NA),
                 panel.background = element_rect(fill = "white", color = NA)),
         
         plot_grid(
           
           video_coords_count |>
             ggplot(aes(longitude, y=latitude)) +
             geom_tile(aes(fill = breaks)) +
             scale_fill_brewer(palette = "YlOrBr", direction = 1) +
             theme_void() +
             labs(
               title = "Geolocated videos 2020-24 (n = 3,398,143)",
               fill = NULL,
             ) +
             guides(fill = guide_legend(nrow = 2)) +
             theme(legend.position = "bottom",
                   plot.background = element_rect(fill = "white", color = NA),
                   panel.background = element_rect(fill = "white", color = NA)),
           
           comments_coords_count |>
             ggplot(aes(longitude, y=latitude)) +
             geom_tile(aes(fill = breaks)) +
             scale_fill_brewer(palette = "YlGnBu", direction = 1) + 
             theme_void() +
             labs(
               title = "Geolocated comments 2020-24 (n = 43,900,491)",
               fill = NULL,
             ) +
             guides(fill = guide_legend(nrow = 2)) +
             theme(legend.position = "bottom",
                   plot.background = element_rect(fill = "white", color = NA),
                   panel.background = element_rect(fill = "white", color = NA)),
           
           nrow = 1
         ),
         
         nrow = 2, rel_heights = c(1, .5)
         
       )

ggsave(filename = "figures/fig-4.jpg", width = 11, height =12,
       fig_4
)

ggsave(filename = "figures/fig-4.pdf", width = 11, height =12,
       fig_4
)
```

## Figure 5

```{r}
library(tidyverse)

load("data/fig-5.RData")

base <- 
  ggplot(daily_range, aes(x = date)) +
  geom_rect(aes(xmin = as.Date("2020-01-01"), xmax = as.Date("2024-06-01"), 
                ymin = period_range$quant25[1], 
                ymax = period_range$quant75[1]), 
            fill = "lightblue", alpha = 0.5) +
  geom_line(aes(y = mean_score), size = .4, alpha = .5) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  scale_y_log10(
    breaks = scales::trans_breaks("log10", function(x) 10^x),
    labels = label_comma() 
  ) +
  scale_x_date(
    date_breaks = "6 months",
    date_labels = "%b %Y"
  ) +
  # Horizontal reference lines
  geom_hline(yintercept = period_range$median[1]) +
  geom_hline(yintercept = 0.05, alpha = .3, linetype = "dashed") +
  # Vertical event lines
  geom_vline(xintercept = as.Date("2020-03-10"), linetype = "dashed", color = "darkgrey", alpha = 0.7) +
  geom_vline(xintercept = as.Date("2020-06-05"), linetype = "dashed", color = "darkgrey", alpha = 0.7) +
  geom_vline(xintercept = as.Date("2020-07-15"), linetype = "dashed", color = "darkgrey", alpha = 0.7) +
  geom_vline(xintercept = as.Date("2023-01-05"), linetype = "dashed", color = "darkgrey", alpha = 0.7) +
  geom_vline(xintercept = as.Date("2024-06-02"), linetype = "dashed", color = "red", alpha = 0.7) +
  # Event annotations with segments pointing to dates
  # March 2020 - Women's protests against violence (MOVED ABOVE SPIKE)
  annotate("segment", 
           x = as.Date("2020-03-10"), xend = as.Date("2020-03-10"),
           y = 0.30, yend = 0.24, 
           arrow = arrow(length = unit(0.15, "cm")), color = "darkgrey") +
  annotate("text", 
           x = as.Date("2020-03-10"), y = 0.35,
           label = "Mar 2020\nMassive protests\nagainst violence", 
           size = 2.5, hjust = 0.5) +
  # June 2020 - Guadalajara protests
  annotate("segment", 
           x = as.Date("2020-06-05"), xend = as.Date("2020-06-05"),
           y = 0.125, yend = 0.055, 
           arrow = arrow(length = unit(0.15, "cm")), color = "darkgrey") +
  annotate("text", 
           x = as.Date("2020-06-05"), y = 0.145,
           label = "Jun 2020\nGuadalajara\nprotests", 
           size = 2.5, hjust = 0.5) +
  # July 2020 - Guanajuato crime wave
  annotate("segment", 
           x = as.Date("2020-07-15"), xend = as.Date("2020-07-15"),
           y = 0.065, yend = 0.055, 
           arrow = arrow(length = unit(0.15, "cm")), color = "darkgrey") +
  annotate("text", 
           x = as.Date("2020-07-15"), y = 0.075,
           label = "Jul 2020\nGuanajuato\ncrime wave", 
           size = 2.5, hjust = 0.5) +
  # January 2023 - Sinaloa unrest
  annotate("segment", 
           x = as.Date("2023-01-05"), xend = as.Date("2023-01-05"),
           y = 0.065, yend = 0.055, 
           arrow = arrow(length = unit(0.15, "cm")), color = "darkgrey") +
  annotate("text", 
           x = as.Date("2023-01-05"), y = 0.08,
           label = "Jan 2023\nSinaloa unrest\n(Guzmán arrest)", 
           size = 2.5, hjust = 0.5) +
  # 2024 Election campaign box
  annotate("rect", 
           ymin = 0.016, ymax = 0.032, 
           xmin = as.Date("2024-01-01"), xmax = as.Date("2024-06-01"),
           fill = "red", alpha = 0.2, color = "red") +
  annotate("text", 
           x = as.Date("2024-03-15"), y = 0.07,
           label = "2024 Federal election\ncampaign (zoomed below)", 
           size = 2.5, hjust = 0.5) +
  labs(
    title = "Average Daily Comment Violence Scores",
    x = NULL, 
    y = NULL
  )
  


base2 <- 
  ggplot(daily_range %>%
           dplyr::filter(date > as.Date("2024-01-01") &
                           date < as.Date("2024-06-03")), aes(x = date)) +
  geom_rect(aes(xmin = as.Date("2024-01-01"), xmax = as.Date("2024-06-01"), 
                ymin = period_range$quant25[1], 
                ymax = period_range$quant75[1]), 
            fill = "lightblue", alpha = 0.5) +
  geom_line(aes(y = mean_score), size = .4, alpha = .5) +
  labs(
    title = "2024 Federal election campaign",
    y = NULL, y = NULL,
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  # Vertical lines
  geom_hline(yintercept = period_range$median[1]) +
  geom_vline(xintercept = as.Date("2024-06-02"), linetype = "dashed", color = "red", alpha = 0.7) +
  annotate("text", 
           x = as.Date("2024-05-28"),  # Center of the box
           y = 0.032,  # Just above the box (ymax = 0.1, so 0.12)
           label = "Election\nday") +
  labs(x = NULL, y = NULL)



ggsave(filename = "figures/fig-5.jpg", width = 12, height = 10,
       
       cowplot::plot_grid(
         
       base +
         scale_x_date(
           limits = c(as.Date("2020-01-01"), as.Date("2024-06-01")),
           date_breaks = "6 months",
           date_labels = "%Y-%m"),
       
       base2 +
         scale_x_date(
           limits = c(as.Date("2024-01-01"), as.Date("2024-06-01")),
           date_breaks = "1 month",
           date_labels = "%Y-%m"),
       
       nrow = 2, rel_heights = c(1,1.2)
       
       )
)

ggsave(filename = "figures/fig-5.pdf", width = 12, height = 10,
       
       cowplot::plot_grid(
         
         base +
           scale_x_date(
             limits = c(as.Date("2020-01-01"), as.Date("2024-06-01")),
             date_breaks = "6 months",
             date_labels = "%Y-%m"),
         
         base2 +
           scale_x_date(
             limits = c(as.Date("2024-01-01"), as.Date("2024-06-01")),
             date_breaks = "1 month",
             date_labels = "%Y-%m"),
         
         nrow = 2, rel_heights = c(1,1.2)
         
       )
)

```

## Figure 6

```{r}
library(tidyverse)
library(readr)
library(sf)
library(cowplot)

priogrid.sf <- 
  read_sf("data/priogrid-shape/priogrid_mex.shp")
priogrid.df <- 
  st_drop_geometry(priogrid.sf)


vpi_acled_homicides_priogrid <- 
  read_csv("data/vpi_acled_homicides_priogrid.csv",
           show_col_types = FALSE)

vpi_acled_homicides_priogrid <- 
  vpi_acled_homicides_priogrid %>%
  dplyr::left_join(priogrid.df, by = "gid")

vpi_acled_homicides_priogrid <- 
  vpi_acled_homicides_priogrid %>%
  dplyr::mutate(vpi_cut = ntile(var1pred_idw_pop_w_std, 9),
                acled_cut = ntile(AC_total_fatalities_pop_w, 9),
                homicide_cut = ntile(crime_count_pop_w, 9),
                vpi_label = factor(vpi_cut, 
                                   levels = c(1:9),
                                   labels = 
                                     c("Q1: Min", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9: Max")),
                acled_label = factor(acled_cut, 
                                   levels = c(1:9),
                                   labels = 
                                     c("Q1: Min", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9: Max")),
                homicide_label = factor(homicide_cut, 
                                   levels = c(1:9),
                                   labels = 
                                     c("Q1: Min", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9: Max")))

fig_6 <- 
  plot_grid(
    
    ggplot(vpi_acled_homicides_priogrid, aes(x = xcoord, ycoord)) +
      geom_tile(aes(fill = vpi_label), size = .2, alpha = .5) +
      scale_fill_brewer(palette = "YlOrRd", direction = 1,
                        na.translate = FALSE) +
      theme_void() +
      labs(
        title = "a. VPI",
        fill = NULL,
      )  +
      guides(fill = guide_legend(nrow = 1)) +
      theme(legend.position = "bottom",
            legend.text = element_text(size = 6)),
    
    ggplot(vpi_acled_homicides_priogrid, aes(x = xcoord, ycoord)) +
      geom_tile(aes(fill = acled_label), size = .2, alpha = .5) +
      scale_fill_brewer(palette = "YlOrRd", direction = 1,
                        na.translate = FALSE) +
      theme_void() +
      labs(
        title = "b. ACLED",
        fill = NULL,
      ) +
      guides(fill = guide_legend(nrow = 1)) +
      theme(legend.position = "bottom",
            legend.text = element_text(size = 6)),
    
    ggplot(vpi_acled_homicides_priogrid, aes(x = xcoord, ycoord)) +
      geom_tile(aes(fill = homicide_label), size = .2, alpha = .5) +
      scale_fill_brewer(palette = "YlOrRd", direction = 1,
                        na.translate = FALSE) +
      theme_void() +
      labs(
        title = "C. Homicides",
             fill = NULL,
           ) +
           guides(fill = guide_legend(nrow = 1)) +
           theme(legend.position = "bottom",
                 legend.text = element_text(size = 6)),
  
  nrow = 1)

ggsave(file = "figures/fig-6.png", width = 15, height = 4,
       fig_6)

ggsave(file = "figures/fig-6.pdf", width = 15, height = 4,
       fig_6)

```

## Figure 7

```{r}
library(tidyverse)
library(readr)
library(cowplot)

acled_month_priogrid <- 
  read_csv("data/acled_month_priogrid.csv", show_col_types = FALSE)

homicide_month_priogrid <- 
  read_csv("data/homicide_month_priogrid.csv")

fig_7 <- 
  plot_grid(
    
    ggplot(acled_month_priogrid, 
           aes(x = AC_total_fatalities, y = `_reghdfe_resid`)) +
      geom_jitter(alpha = .5) + 
      labs(x = "Fatalities (iym)", y = "Residuals (iym)",
           title = "a. ACLED"), 
    
    
    ggplot(homicide_month_priogrid, 
           aes(x = crime_count, y = `_reghdfe_resid`)) +
      geom_jitter(alpha = .5) + 
      labs(x = "Fatalities (iym)", y = "Residuals (iym)",
           title = "b. Homicides"), 
    
    ggplot(acled_month_priogrid, 
           aes(x = marg_index_2020, y = `_reghdfe_resid`)) +
      geom_jitter(alpha = .5) + 
      labs(x = "Marginalization (2020)", y = "Residuals (iym)",
           title = "c. ACLED"), 
    
    ggplot(homicide_month_priogrid, 
           aes(x = marg_index_2020, y = `_reghdfe_resid`)) +
      geom_jitter(alpha = .5) + 
      labs(x = "Marginalization (2020)", y = "Residuals (iym)",
           title = "d. Homicides"), 
    
    
    nrow = 2, ncol = 2)

ggsave(file = "figures/fig-7.jpg", width = 12, height = 12,
       fig_7)

ggsave(file = "figures/fig-7.pdf", width = 12, height = 12,
       fig_7)

```

## Figure 8

```{r}
library(readr)
library(sf)
library(tidyverse)

priogrid.sf <- 
  read_sf("data/priogrid-shape/priogrid_mex.shp")
priogrid.df <- 
  st_drop_geometry(priogrid.sf)

acled_resid_priodgrid <- 
  read_csv("data/acled_resid_priodgrid.csv", show_col_types = FALSE)

homicide_resid_priodgrid <- 
  read_csv("data/homicide_resid_priodgrid.csv", show_col_types = FALSE)

marg_index_2020_priogrid <- 
  read_csv("data/marg_index_2020_priogrid.csv", show_col_types = FALSE)

this_dat <- 
  acled_resid_priodgrid %>%
  dplyr::left_join(homicide_resid_priodgrid, by = "gid") %>%
  dplyr::left_join(marg_index_2020_priogrid, by = "gid") %>%
  dplyr::left_join(priogrid.df, by = "gid")

this_dat <-
  this_dat %>%
  dplyr::mutate(marg_cut = ntile(IM_2020, 5),
                acled_cut = ntile(ACLED_reghdfe_resid_abs, 5),
                homicide_cut = ntile(crime_reghdfe_resid_abs, 5),
                marg_label = factor(marg_cut, 
                                   levels = c(1:5),
                                   labels = 
                                     c("Q1: Min", "Q2", "Q3", "Q4", "Q5: Max")),
                acled_label = factor(acled_cut, 
                                   levels = c(1:5),
                                   labels = 
                                     c("Q1: Min", "Q2", "Q3", "Q4", "Q5: Max")),
                homicide_label = factor(homicide_cut, 
                                   levels = c(1:5),
                                   labels = 
                                     c("Q1: Min", "Q2", "Q3", "Q4", "Q5: Max")))


fig_8 <-
  plot_grid(
    
    ggplot(this_dat, aes(x = xcoord, y = ycoord)) +
      geom_tile(aes(fill = acled_label)) +
      scale_fill_brewer(palette = "YlOrRd", direction = 1,
                        na.translate = FALSE) +
      theme_void() +
      labs(
        title = "a. ACLED (error term)",
        fill = NULL,
      ) +
      guides(fill = guide_legend(nrow = 1)) +
      theme(legend.position = "bottom"),
    
    ggplot(this_dat, aes(x = xcoord, y = ycoord)) +
      geom_tile(aes(fill = homicide_label)) +
      scale_fill_brewer(palette = "YlOrRd", direction = 1,
                        na.translate = FALSE) +
      theme_void() +
      labs(
        title = "b. Homicide (error term)",
        fill = NULL,
      ) +
      guides(fill = guide_legend(nrow = 1)) +
      theme(legend.position = "bottom"),
    
    ggplot(this_dat, aes(x = xcoord, y = ycoord)) +
      geom_tile(aes(fill = marg_label)) +
      scale_fill_brewer(palette = "YlOrRd", direction = 1,
                        na.translate = FALSE) +
      theme_void() +
      labs(
        title = "c. Marginalization (index)",
        fill = NULL,
      ) +
      guides(fill = guide_legend(nrow = 1)) +
      theme(legend.position = "bottom"),
    
    
    ncol = 3
    
  )

ggsave(file = "figures/fig-8.png", width = 15, height = 4,
       fig_8)

ggsave(file = "figures/fig-8.pdf", width = 15, height = 4,
       fig_8)
  
```
